Protocol 2 - Petrographic data

The following example applies protocol 2 to confirm workshops’ petrographic groups.

Protocol 2 consist in:

  1. Select ordinal petrographic data;
  2. Transform to ranks;
  3. Extended Gower distance, using:
  1. Relative ranking difference (RRD)
  2. Neighbor interchange (NI)
  1. Apply ordination procedure:
  1. Principal Coordinates Analysis (PCoA)
  2. Non-metric Dimensional Scaling (NMDS)
  1. Perform PERMANOVA & PERMDISP tests;

Last, search for outliers and re-do protocol excluding outliers.

NOTE: The initial procedures must be ran at least once before any protocol can be applied.

The key references on the Extended Gower distance are:

Pavoine, S., Vallet, J., Dufour, A.-B., Gachet, S., Daniel, H., 2009. On the challenge of treating various types of variables: application for improving the measurement of functional diversity. Oikos 118, 391-402. doi:10.1111/j.1600-0706.2008.16668.x

Podani, J., 1999. Extending Gower’s General Coefficient of Similarity to Ordinal Characters on JSTOR. Taxon 48, 331-340. doi:10.2307/1224438

Gower, J.C., 1971. A General Coefficient of Similarity and Some of Its Properties. Biometrics 27, 857-871. doi:10.2307/2528823

Ordination procedure

Depending on which type of distance calculation (RRD/NI), protocol 2 performs different ordination methods (PCoA/NMDS). Both PCoA and NMDS require specifying the number of dimensions in which to project the data. Therefore, you must generate specific 2D and 3D ordination objects:

prot2a_2d <- apply_ordination(cleanAmphorae[!isShipwreck,],
                              "2a", # select protocol 2a (RRD & PCoA)
                              exception_columns = excep_cols,
                              variable_tags = varCode)

prot2b_2d <- apply_ordination(cleanAmphorae[!isShipwreck,],
                              "2b", # select protocol 2a (NI & NMDS)
                              exception_columns = excep_cols,
                              variable_tags = varCode)

prot2a_3d <- apply_ordination(cleanAmphorae[!isShipwreck,],
                              "2a", # select protocol 2a (RRD & PCoA)
                              exception_columns = excep_cols,
                              variable_tags = varCode,
                              dimensions = 3)

prot2b_3d <- apply_ordination(cleanAmphorae[!isShipwreck,],
                              "2b", # select protocol 2a (NI & NMDS)
                              exception_columns = excep_cols,
                              variable_tags = varCode,
                              dimensions = 3)

The ordination objects generated with protocol 2 are different from those in protocol 1 since it uses different functions. However, the main components are still the same: the projection of observations or scores (points) and of variables or loadings.

class(prot2a_2d)
#> [1] "list"
names(prot2a_2d)
#>  [1] "points"        "eig"           "x"             "ac"           
#>  [5] "GOF"           "sub2D"         "GOF2_2D"       "loadings"     
#>  [9] "variable_tags" "name"          "dist_matrix"
class(prot2b_2d)
#> [1] "metaMDS" "monoMDS"
names(prot2b_2d)
#>  [1] "nobj"          "nfix"          "ndim"          "ndis"         
#>  [5] "ngrp"          "diss"          "iidx"          "jidx"         
#>  [9] "xinit"         "istart"        "isform"        "ities"        
#> [13] "iregn"         "iscal"         "maxits"        "sratmx"       
#> [17] "strmin"        "sfgrmn"        "dist"          "dhat"         
#> [21] "points"        "stress"        "grstress"      "iters"        
#> [25] "icause"        "call"          "model"         "distmethod"   
#> [29] "distcall"      "distance"      "converged"     "tries"        
#> [33] "engine"        "species"       "data"          "init_seed"    
#> [37] "trymax"        "sub_stress"    "sub2D"         "GOF2_2D"      
#> [41] "loadings"      "variable_tags" "name"          "dist_matrix"

Test the given fabric groups

The fabric groups defined in previous studies can be tested against Protocol 2 distance matrices. We can use either “prot2a_2d$dist_matrix” or “prot2a_3d$dist_matrix”, because they are the same. Remember that this test batch may take several minutes.

prot2a_tests <- test_groups(prot2a_2d$dist_matrix, 
                            factor_list$FabricGroup)
#> [1] "initiating test batch..."
#> [1] "vegan::anosim done."
#> [1] "vegan::betadisper done."
#> [1] "vegan::permutest done."
#> [1] "vegan::adonis done."
#> [1] "Test batch completed."
prot2b_tests <- test_groups(prot2b_2d$dist_matrix,
                            factor_list$FabricGroup)
#> [1] "initiating test batch..."
#> [1] "vegan::anosim done."
#> [1] "vegan::betadisper done."
#> [1] "vegan::permutest done."
#> [1] "vegan::adonis done."
#> [1] "Test batch completed."

These tests were explained in protocol 1.

Biplots

The details on how to create biplots is already explained in protocol 1. Concerning protocol 2, we can compare the results of version 2a (RRD, PCoA) and 2b (NI, NMDS).

Biplot 2D

arrows_label_adj <- rbind(c(.5,.8),c(.5,1),c(.5,1),c(.5,0),c(.5,1),
                          c(.5,0),c(0,.5))
row.names(arrows_label_adj) <- c("L48","L24","L5","L36","S7",
                                 "S8","S11")

biplot2d3d::biplot_2d(prot2a_2d,
                      ordination_method = "PCoA",
                      invert_coordinates = c (TRUE,TRUE),
                      xlim = c(-.26,.35),
                      ylim = c(-.31,.35),
                      point_type = "point",
                      groups = factor_list$FabricGroup,
                      group_color = color_list$FabricGroup,
                      group_label_cex = 0.6,
                      arrow_mim_dist = 0.5,
                      arrow_label_cex = 0.6,
                      arrow_fig = c(.6,.95,0,.35),
                      arrow_label_adj_override = arrows_label_adj,
                      subtitle = prot2a_2d$sub2D,
                      test_text = prot2a_tests$text(prot2a_tests),
                      test_cex = 0.8,
                      test_fig = c(0, 0.5, 0.62, .99),
                      fitAnalysis_fig = c(0,.7,.05,.5),
                      output_type = "preview")
protocol 2a

protocol 2a

arrows_label_adj <- rbind(c(.5,1),c(.5,0),c(.5,1),c(.5,1),c(.5,0),
                          c(0,.5),c(1,.5))
row.names(arrows_label_adj) <- c("S7","S8","CLAY","L24","L43",
                                 "L5","L36")

biplot2d3d::biplot_2d(prot2b_2d,
                      ordination_method = "NMDS",
                      xlim = c(-.42,.38),
                      ylim = c(-.45,.25),
                      point_type = "point",
                      groups = factor_list$FabricGroup,
                      group_color = color_list$FabricGroup,
                      group_label_cex = 0.6,
                      arrow_mim_dist = .5,
                      arrow_label_cex = 0.6,
                      arrow_fig = c(.6,.95,0,.35),
                      arrow_label_adj_override = arrows_label_adj,
                      subtitle = prot2b_2d$sub2D,
                      test_text = prot2b_tests$text(prot2b_tests),
                      test_cex = 0.8,
                      test_fig = c(0, 0.5, 0.62, .99),
                      fitAnalysis_stress_axis_cex = 0.8,
                      fitAnalysis_fig = c(.1,.6,.1,.4),
                      output_type = "preview")
protocol 2b

protocol 2b

Biplot 3D

biplot2d3d::biplot_3d(prot2a_3d,
                       ordination_method = "PCoA",
                       point_type = "point",
                       groups = factor_list$FabricGroup,
                       group_color = color_list$FabricGroup,
                       group_representation = "stars",
                       star_centroid_radius = 0,
                       star_label_cex = .8,
                       arrow_min_dist = .5,
                       arrow_body_length = .025,
                       subtitle = prot2a_3d$sub3D,
                       test_text = prot2a_tests$text(prot2a_tests),
                       test_cex = 1.25,
                       test_fig = c(0, 0.5, 0.65, .99),
                       view_zoom = 0.9)

biplot2d3d::animation(directory = directories$prot2,
                       file_name = "Prot2a_Biplot3D")
Prot2a_Biplot3D.gif

Prot2a_Biplot3D.gif

biplot2d3d::biplot_3d(prot2b_3d,
                       ordination_method = "NMDS",
                       point_type = "point",
                       groups = factor_list$FabricGroup,
                       group_color = color_list$FabricGroup,
                       group_representation = "stars",
                       star_centroid_radius = 0,
                       star_label_cex = .8,
                       arrow_min_dist = .5,
                       arrow_body_length = .025,
                       subtitle = prot2b_3d$sub3D,
                       test_text = prot2b_tests$text(prot2b_tests),
                       test_cex = 1.25,
                       test_fig = c(0, 0.5, 0.65, .99),
                       view_zoom = 0.9)

biplot2d3d::animation(directory = directories$prot2,
                       file_name = "Prot2b_Biplot3D")
Prot2b_Biplot3D_snapshot.png

Prot2b_Biplot3D_snapshot.png

Back to index